restart;
with(LinearAlgebra):


#3D unsteady velocity field U3 in first-order Taylor approximation
uuf := uu + uu_x*x + uu_y*y + uu_z*z + uu_t*t;
vvf := vv + vv_x*x + vv_y*y + vv_z*z + vv_t*t;
wwf := ww + ww_x*x + ww_y*y + ww_z*z + ww_t*t;
U3 := Vector([uuf,vvf,wwf]);


#current particle velocity
V3 := Vector([u,v,w]);

#gravity Vector
G3 := Vector([ug,vg,wg]);


#3D Spatial Jacobian matrix J3:
J3 := Matrix([
[diff(uuf,x),diff(uuf,y),diff(uuf,z)],
[diff(vvf,x),diff(vvf,y),diff(vvf,z)],
[diff(vvf,x),diff(wwf,y),diff(wwf,z)]]);


#7D vector field P7, defined in (19)
dx := u;
dy := v;
dz := w;
du := (uuf-u)/r + ug;
dv := (vvf-v)/r + vg;
dw := (wwf-w)/r + wg;
dt := 1;
P7 := Vector([dx,dy,dz,du,dv,dw,dt]);


#6D vector field U6, 
dx := u;
dy := v;
dz := w;
du := (uuf-u)/r + ug;
dv := (vvf-v)/r + vg;
dw := (wwf-w)/r + wg;
U6 := Vector([dx,dy,dz,du,dv,dw]);


#7D Jacobian matrix J7 of P7, (20)
J7 := Matrix([
[diff(dx,x),diff(dx,y),diff(dx,z),diff(dx,u),diff(dx,v),diff(dx,w),diff(dx,t)],
[diff(dy,x),diff(dy,y),diff(dy,z),diff(dy,u),diff(dy,v),diff(dy,w),diff(dy,t)],
[diff(dz,x),diff(dz,y),diff(dz,z),diff(dz,u),diff(dz,v),diff(dz,w),diff(dz,t)],
[diff(du,x),diff(du,y),diff(du,z),diff(du,u),diff(du,v),diff(du,w),diff(du,t)],
[diff(dv,x),diff(dv,y),diff(dv,z),diff(dv,u),diff(dv,v),diff(dv,w),diff(dv,t)],
[diff(dw,x),diff(dw,y),diff(dw,z),diff(dw,u),diff(dw,v),diff(dw,w),diff(dw,t)],
[diff(dt,x),diff(dt,y),diff(dt,z),diff(dt,u),diff(dt,v),diff(dt,w),diff(dt,t)]
]);


#6D Jacobian matrix J6 of U6 (9)
J6 := 
Matrix([
[diff(dx,x),diff(dx,y),diff(dx,z),diff(dx,u),diff(dx,v),diff(dx,w)],
[diff(dy,x),diff(dy,y),diff(dy,z),diff(dy,u),diff(dy,v),diff(dy,w)],
[diff(dz,x),diff(dz,y),diff(dz,z),diff(dz,u),diff(dz,v),diff(dz,w)],
[diff(du,x),diff(du,y),diff(du,z),diff(du,u),diff(du,v),diff(du,w)],
[diff(dv,x),diff(dv,y),diff(dv,z),diff(dv,u),diff(dv,v),diff(dv,w)],
[diff(dw,x),diff(dw,y),diff(dw,z),diff(dw,u),diff(dw,v),diff(dw,w)]
]);


x := 0;
y := 0;
z := 0;
t := 0;


#compare (20) in the paper
J7;

#the field in (7) in the paper
F3 := Vector([ufff,vfff,wfff]);
F6 := Vector([ufff,vfff,wfff,0,0,0]);


#compare (11) in the paper:
W := (U3-V3)/r + G3;

#compare (25) in the paper:
Multiply( J6, U6-F6 );


#search Eigenvector E7 of J7 belonging to Eigenvalue 0:
E7 := Vector([ufff,vfff,wfff,0,0,0,1]);

sol1 := solve({
Multiply(J7,E7)[1]=0,
Multiply(J7,E7)[2]=0,
Multiply(J7,E7)[3]=0,
Multiply(J7,E7)[4]=0,
Multiply(J7,E7)[5]=0,
Multiply(J7,E7)[6]=0,
Multiply(J7,E7)[7]=0},{ufff,vfff,wfff});
assign(sol1);

#This is exactly the vector in (7):
F3 := Vector([ufff,vfff,wfff]);



#compare (9) in the paper:
J6;

#compare (11) in the paper:
W := (U3-V3)/r + G3;

#compare (12) in the paper:
Multiply(J6,U6);